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We present an algorithm for solving the self-consistency equations of the dynamical mean-field 
theory (DMFT) with high precision and efficiency at low temperatures. In each DMFT iteration, 
the impurity problem is mapped to an auxiliary Hamiltonian, for which the Green function is 
computed by combining determinantal quantum Monte Carlo (BSS-QMC) calculations with a 
multigrid extrapolation procedure. The method is numerically exact, i.e., yields results which are 
free of significant Trotter errors, but retains the BSS advantage, compared to direct QMC impurity 
solvers, of linear (instead of cubic) scaling with the inverse temperature. The new algorithm is 
applied to the half-filled Hubbard model close to the Mott transition; detailed comparisons with 
exact diagonalization, Hirsch-Fye QMC, and continuous-time QMC are provided. 

PACS numbers: 02.70.Ss, 71.10.Fd, 71.27.-|-a 



I. INTRODUCTION 

The dynamical mean- field theory (DMFT)[ll-I3 and 
its cluster extensions Q are powerful approaches for 
the numerical treatment of correlated electron systems, 
both in the model context and for materials science, e., 
embedded in the LDA-fDMFT or GW-hDMFT ' 
[llj frameworks which extend density functional theory 
to strongly correlated materials [1, Il2j . Recently, many 
DMFT studies have also appeared in the context of 
ultracold fcrmions on optical lattices [l3l - [l5| . The 
DMFT reduces electronic lattice models to impurity 
problems, which have to be solved self-consistently |16l - 
Il8| . A challenging part of this iterative procedure is 
the computation of the interacting Green function for 
a given impurity configuration (defined by the fixed local 
interactions and the self-consistent Weiss field). Thus, 
the availability of efficient and reliable impurity solvers 
determines the complexity of models and the parameter 
space that can be accessed using the DMFT. 

Quantum Monte Carlo (QMC) impurity solvers al- 
low for numerically exact solutions of the DMFT self- 
consistency equations at finite temperatures. In the 
case of the Hirsch-Fye auxiliary field (HF-QMC) method 
[igI [igl . [20| , all raw estimates contain systematic errors 
due to the inherent Trotter decomposition and associated 
imaginary-time discretization [l^ [2l|; unbiased results 
can only be obtained after an extrapolation of the dis- 
cretization interval Ar — >■ [lHjli^l- Diagrammatic QMC 
impurity solvers (23 - [27j sample partition function and 
Green functions in continuous (imaginary) time (CT), 
i.e., avoid systematic biases. However, in all of these 
direct QMC approaches, the computational effort scales 
cubically [28| with the inverse temperature /3 ~ l/fcaT', 
which limits their access to low-temperature phases. 

Exact diagonalization (ED) based impurity solvers [2^ 
require a discrete representation of the impurity action in 
terms of an auxiliary Hamiltonian, which is then solved 



either by full diagonalization (for evaluations at arbitrary 
temperature) or using a Lanczos procedure [s^l (e.g., at 
T = 0). As the numerical effort scales exponentially with 
the number Nf, of auxiliary "bath" sites, Ni, has to be 
kept quite small, which introduces, again, a bias and is a 
particularly severe limitation for multi-orbital or cluster 
DMFT studies at finite temperatures. 

Recently, Khatami et al. proposed another 
Hamiltonian-based scheme [sif . in which the Green 
function and other relevant properties of the auxiliary 
problem are computed using the determinantal BSS- 
QMC method developed by Blankenbecler, Scalapino, 
and Sugar [s^l- The advantage of this scheme, compared 
to ED, is the possibility of using more bath sites (due 
to cubic instead of exponential scaling with Nb); the 
advantage over the direct QMC impurity solvers is the 
linear, instead of cubic, scaling in /3 (SJl- The authors 
established the feasibility of the method and proved 
that the associated sign problem (arising at general 
band filling in cluster extensions of DMFT, in frustrated 
lattices, and for generic multiband models) converges to 
that of HF-QMC for sufficiently fine bath discretization 
|3l| . However, as all BSS-QMC applications to date, the 
Green functions and all observable estimates resulting 
from their implementation suffer from systematic Trotter 
errors. 

In this work, we construct a similar algorithm where 
the Trotter bias inherent in the BSS Green functions is 
eliminated using a multigrid procedure before feeding 
them back in the self-consistency cycle. As a DMFT 
building block, the resulting method is an exact quasi- 
CT QMC impurity solver with linear scaling in the 
inverse temperature. Its scaling advantage over direct 
QMC impurity solvers should allow access to lower 
temperatures, in particular in multi-orbital and cluster 
DMFT studies. 

The paper is organized as follows: In Sec. |TI1 we 
briefly review the DMFT equations and the BSS-QMC 
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algorithm and fully specify our multigrid DMFT-BSS 
approach. As a test case, the new method is applied 
(in single-site DMFT) to the half-filled Hubbard model 
in the vicinity of the Mott transition in Sec. IIIII Here, 
we first focus on the Green function at moderately 
low temperature T = t/25 and then discuss important 
obscrvablcs, namely the double occupancy and quasipar- 
ticle weight, also at lower temperatures. The accuracy 
of our approach is established by comparisons with the 
results of the (multigrid) HF-QMC, ED and CT-QMC 
solvers, as well as with the previous (finite Ar) DMFT- 
BSS implementation. We show that our elimination of 
the Trotter error improves the results dramatically. We 
also discuss the impact of the bath discretization and 
establish convergence to the thermodynamic limit. A 
summary and outlook conclude the paper in Sec. IIVI 



II. THEORY AND ALGORITHMS 

In this section, we lay out the proposed algorithm 
for solving the DMFT self-consistency equations without 
significant Trotter errors and with a computational effort 
that grows only linearly with the inverse temperature. 
We start out by reviewing the general DMFT framework 
and established methods (ED, HF-QMC, CT-QMC) for 
its solution in Sec. IH Al in sufficient detail to expose 
the similarities and differences with respect to the new 
method. Here we also discuss some algorithmic choices, 
in particular regarding the Hamiltonian representation 
in our ED implementation, that are essential ingredients 
also for the BSS-QMC based approaches. We then turn 
to the BSS-QMC method and its applicability in the 
DMFT context in Sec. IH Bl and specify, finally, our 
new numerically exact implementation in Sec. IH CI For 
simplicity, and in line with the numerical results to be 
presented in Sec. IHII we write down the formalism for 
the single-band Hubbard model and the original, single- 
site variant of DMFT. Extensions to cluster DMFT (and 
to multiband models) should be straightforward, but 
require some generalizations (e.g. for the treatment of 
offdiagonal Green functions) and will be pursued in a 
subsequent publication. 



A. DMFT and established impurity solvers 

The Hubbard model on a lattice or graph - We consider 
the single-band Hubbard model 

H = Ho + Fint = t,j c\^c^^ + rij^rijj, , (1) 



where {Cia) creates (annihilates) an electron with 
spin a € {ti 4-} on lattice site i; ni„ = c\^c^^ is the 
corresponding density, tij = tji the hopping amplitude 
between sites i and j (or the local potential for i = j); U 
quantifies the on-site interaction. Usually, the hopping 




FIG. 1. (Color online) Mapping of the original lattice problem 
(a), with local interaction U and hopping t, on a single impu- 
rity (b), embedded in an effective bath Q. (c) Discretization 
of the bath in terms of an auxiliary Hamiltonian (treatable 
with ED or BSS-QMC), here with star topology. 

is defined to be translationally invariant, e.g., ty = —t 
for nearest-neighbor bonds on an infinite mathematical 
lattice [as illustrated in Fig. [TJa) for a square lattice] 
or on a finite cluster with periodic boundary conditions. 
However, neither the DMFT nor direct QMC approaches 
to the Hubbard model depend crucially on such assump- 
tions, as will be discussed in Sec. IHBI 

General DMFT self-consistency procedure - If all 
lattice sites are equivalent and for spatially homogeneous 
phases, the DMFT maps the original lattice problem ([T|), 
illustrated in Fig. [TJa), onto a single- impurity Anderson 
model [Fig.lUb)], which has to be solved self-consistently. 
The impurity problem is defined by its action 
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(2) 



" 
-U J dT^^T)lP^{T)lPl{T)^^{T) , 



here in imaginary time t € [0,/3] and in terms of 
Grassmann fields ip, ip* ■ Q is the "bath" Green function, 
i.e., the noninteracting Green function of the impurity, 
which is related to the full impurity Green function G, 



G.(r) = -(r.^,(r)V:(0)). 



(3) 



(with the time ordering operator %), and the self-energy 
E by the (impurity) Dyson equation 



G^^(iw„) = g„'^{iujn) - 'E^^iujn) 



(4) 



here written in terms of fermionic Matsubara frequencies 
a;„ — {2n + l)TTT at finite temperature T; here and in the 
following, we set ft = fee = 1- 

The central DMFT assumption is that of a local self- 
energy on the lattice T,ij„{iujn) = SijT,ua{iuJn), which 
is identified with the impurity self-energy. Similarly, 
the impurity Green function is identified with the local 
component of the lattice Green function; 

Ga{iuJn) = Gua{iLOn) = {l + [ii^n + fJ- - T:„{iuJn)] 



de 



ioJn + IJ,- T.a(iUJn) 



(5) 



3 



(a) 



A 

4 



lattice Dyson eq. (5) 



impurity Dyson cq. (4) 



impurity problem (3), (2) 



Z^G = G [S; U, T]^<: 




A (b) 



HF-QMC 




FIG. 2. (Color online) (a) Scheme of the general DMFT self-consistency cycle, including the "impurity problem" (dashed box). 
Established impurity solvers include (b) the Hirsch-Fye (HF-QMC) algorithm and (c) exact diagonalization (ED): cf. the main 
text, (d) The proposed algorithm approximates the bath Green function Q in terms of the parameters Vi, of an auxiliary 
Hamiltonian ^ with Ni, "bath" sites [like ED (c)]. Corresponding Green functions are computed using BSS-QMC for a grid 



< Ar, < Arn 



of Trotter discretizations. The subsequent extrapolation of Ar — > yields the Green function, free of 



significant Trotter errors and continuous in r, which is easily Fourier transformed and fed back into the self-consistency cycle. 



where the last expression is valid in the homogeneous 
case, p{e) denotes the noninteracting density of states, 
and It is the matrix with elements tij. 

The general DMFT iteration scheme is illustrated in 
Fig. [2la): starting, e.g., with an initial guess E = Eq 
of the self-energy, the Green function G is computed 
using the lattice Dyson equation ([5]). In a second step, E 
and G yield the bath Green function Q via the impurity 
Dyson equation which defines, in combination with 
the local interactions, the impurity problem [illustrated 
in Fig. [IJb)], the solution of which is the nontrivial part 
of the algorithm. A second application of the impurity 
Dyson equation (jl]), to the resulting G and to Q, yields 
a new estimate of the self-energy E, which closes the 
self-consistency cycle. In the following, we discuss the 
primary options for addressing the impurity problem. 

Direct impurity solvers - One class of methods directly 
evaluates the path integral representation of the Green 
function [Eqs. ([3]) and ([2|)] for a continuous bath Q, 
which corresponds to a DMFT solution of the original 
lattice problem in the thermodynamic limit after self- 
consistency. We will refer to such methods as "direct 
impurity solvers." 

For a long time, the Hirsch-Fye QMC (HF-QMC) 
algorithm has been the method of choice for nonper- 
turbative DMFT calculations HF-QMC is based 

on a discretization of the imaginary time r S [0, /3] 
into A "time shoes" of width Ar = /3/A, a Trotter 
decomposition of the interaction and kinetic terms in Eq. 

and a Hubbard- Stratonovich transformation, which 
replaces the electron-electron interaction by an auxiliary 
binary field on each time slice; the resulting problem is 
then solved employing Wick's theorem and Monte Carlo 
importance sampling over the field configurations. As 
configurations can be updated in the case of a single 
spin flip (i.e., an auxiliary- field change on a single time 
slice) with a matrix- vector operation of cost 0{A^) and 
A local updates are needed for a global configuration 
update, the numerical cost of the HF-QMC algorithm 
scales as A'^. All HF-QMC results have statistical errors 



(which decay as N~^^^ for N "sweeps," each consisting of 
A attempted single-spin updates) and systematic errors 
resulting from the Trotter decomposition. As Ar has to 
be kept constant for roughly constant systematic error 
upon variation of T, the numerical effort of HF-QMC 
scales as the cube of the inverse temperature, (3^. This is 
also true for the numerically exact (unbiased) "multigrid" 
HF-QMC method [H]. 

The integration of the (conventional) HF-QMC 
method into the DMFT self-consistency cycle is illus- 
trated in Fig. [21b) [as a specification of the lower dashed 
box in Fig.[51a)]: a fixed choice of Ar (diamond-shaped 
selection box) defines the grid rj = IAt with < Z < A 
for a Fourier transform (square box) of the Matsubara 
bath Green function Q{iuin) (with |a;„| < Wmax for some 
cutoff frequency Wmax) to the imaginary-time equivalent 
{Q{Ti)}t=o- After apphcation of the HF-QMC algorithm 
(rounded box), the result {G(r/)}^Q is transformed back 
to Matsubara frequencies (square box) ; this step requires 
special care in order to get around the Nyquist theorem, 
e.g. using analytic weak-coupling results ^23, 35, 3^. 

More recently, conceptionally different QMC ap- 
proaches have been formulated, which arc based on 
diagrammatic expansions of the action ^ in continuous 
imaginary time, either in the interaction U (CT-INT (37| ) 
or in the bath hybridization (CT-HYB HI ill), and on 
a stochastic sampling of Feynman dia gram s; CT-AUX 
[Hi is related to the HF-QMC method^. Ah of these 
continuous time (CT-QMC) algorithms require Fourier 
transforms, before and (with the exception of CT-INT) 
after the QMC part; the numerical cost is associated 
primarily with matrix updates, similar to those arising 
in HF-QMC, with a total scaling of the computational 
effort, again, as P'^ . 

Thus, all direct QMC-based impurity solvers are very 
costly at low T, which limits their access to low- 
temperature phases of particular physical interest. 

Auxiliary Hamiltonian and exact diagonalization - 
Another class of numerical approaches, such as the "exact 
diagonalization" methods, cannot directly be applied to 
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the action-based formulation of the impurity problem, 
but requires a Hamiltonian representation [39[. One 
possibility is the "star topology" illustrated in Fig. [UJe) , 
where a central "impurity" site (with the same inter- 
actions as the impurity problem, here U) is connected 
by hopping matrix elements Via to a number Nb of 
nonintcracting "bath sites," each characterized by a local 
potential ei^. In general, this representation has to be 
spin-dependent, leading to the Anderson Hamiltonian 



-ff And = £o + t/n-j. 



+ h.c. 



, (6) 



where (c^) creates (annihilates) an electron with spin 
cr S {t>i} on the impurity site and i^ia) creates 
(annihilates) an electron with spin a on bath site i; 
Tier = c^c^, Tiia = al^a^^ are the corresponding number 
operators. In this work, we consider only nonmagnetic 
phases, which implies spin symmetric bath parameters 
= Vi^, en = ti^. 
For a fixed choice of A^f,, the bath parameters eicr^Via 
are determined such that the noninteracting impurity 
Green function QAnd.a associated with i?And, with 



the Newton method, based on analytic expressions for the 
derivative Vx^ with respect to the bath parameters. Due 
to the multidimensional character of the problem, this 
deterministic method is often trapped in local minima; 
thus, a naive implementation of Newton based methods 
will, in general, not find globally optimal parameters, 
which can induce unphysical fixed points in the DMFT 
iteration procedure. Therefore, we use not only the 
solution {Vi, Ei} of the previous iteration as initialization, 
but perform a large number (up to 1000) of independent 
Newton searches, starting also from random initial pa- 
rameters. Of the resulting locally optimal solutions, we 
choose the one with minimum as the final result of 
the minimization procedure; typically, about 1% of the 
individual searches come close to this (estimated) global 
optimum. 

An advantage of ED, compared to QMC algorithms, 
is that Green functions and spectra can be computed 
directly on the real axis, without analytic continuation; 
however, numerical broadening of the resulting discrete 
peaks is required. This discretization problem is partic- 
ularly severe as the numerical effort of the matrix diago- 
nalization scales exponentially with the total number of 
sites (here A''f, -I- 1), which limits the applicability of ED 
for cluster extensions of DMFT or multiband models. 



?And,, 



(^) 



T/2 



(7) 



is "close" to the target bath Green function Q according 
to some metric (sec below). Note that the resulting 
spectrum ~—lmQxnd,<y{^ + ^0"'') is necessarily discrete 
(for finite A^t,), in contrast to piece- wise smooth spectrum 
of the true bath Green function; in this sense, the 
mapping to a Hamiltonian implies a "bath discretization" 
in frequency space; this step clearly introduces a bias 
which has to be controlled particularly carefully within 
iterative procedures such as the DMFT. 

The integration of this type of approach in the DMFT 
cycle is illustrated for the case of ED in Fig. [2jc): for a 
fixed choice of A";, (diamond shaped box), the parameters 
Vi^ Ci (here and in the following we suppress spin indices) 
are adjusted (rounded box) as to minimize the bath misfit 



X 



= E^" |5And(jw„; {K:,£i}) -~Q{i^^n)Y 



ri=0 



with a cutoff Matsubara frequency iw„^ and the weight- 
ing factor Wn, which can be used to optimize the bath 
parametrization [4^ and which we set to w„ = 1 41 1 . As 



this fit is performed directly on the Matsubara axis, no 
Fourier transform is needed for Q. Using ED (rounded 
box), the Green function G can be evaluated on the 
Matsubara axis [1^; therefore, the DMFT cycle is closed 
without any Fourier transform. 

The minimization of [as defined in Eq. ([S])] is 
performed in our ED and BSS-DMFT calculations using 



Principles of the BSS-QMC algorithm and 
application as a DMFT impurity solver 



In Eq. ([H]), we have used the conventional notation for 
the auxiliary Hamiltonian that emphasizes its interpre- 
tation as an impurity model, e.g., with different creation 
operators for electrons on the central "impurity" site 
(c^) and on the bath sites (a|^), respectively. However, 



and Via ^oi' essentially reproduces the Hubbard 



with the changes 

model (HI) on a graph, just with nonuniform interaction 
{U acting only on site 0) and, possibly, spin-dependent 
hopping amplitudes and local energies. 

As a consequence, the model (O is not only treatable 
with the universal ED approach, but also with more 
specific methods developed for Hubbard-type models. As 
pointed out recently by Khatami et al. [3l| , this includes 
the determinantal quantum Monte Carlo approach by 
Blankenbecler, Scalapino, and Sugar [3^ |42| (BSS- 
QMC), which, thereby, becomes applicable as a DMFT 
impurity solver. In the following, we will first sketch 
the established BSS-QMC approach (for an extended 
discussion, including issues of parallelization, see Ref.l43l) 
and then discuss its application in the DMFT context. 



Similarly to the HF-QMC method (cf. Sec. lEI, 
the BSS-QMC approach is based on a Trotter-Suzuki 
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decomposition, here of the partition function 



Z =Tr 



-I3{Hk+Hv) 



• Zat = Tr 



(9) 
(10) 



where Hy {Hk) corresponds to the interaction (kinetic 
and local potential) contribution to the Hubbard type 
models ([T|) or (|6]) and Ar = /3/A. Again, a discrete 
Hubbard-Stratonovich transformation replaces the inter- 
action term by a binary auxiliary field {/i} with hi{l) = 
±1 at each site i and time slice The trace in Eq. PH)) 
then simplifies to 



^det 

{h} 



M. 



{h}' 



det 



M 



{h} 



with 



(11) 



A/i''> = 1 + Ba.. [{h,{A)}l,] . . . [{h,{l)}li] , 
where B is defined in terms of the hopping matrix K: 
Bi,a [{h^{l)]f=l] - e'^^diag[/n(0. (01 e-^^^. (12) 



The interaction strength is encoded in the parameter A = 
cosh~^(e'^'^'^/^). The computation of thermal averages of 
physical observables O takes the form: 



{h} 

-uih} _ 1 
'"^ At ry 

^At 



(13) 



det 



M. 



{h} 



det 



M 



At particle-hole symmetry, the weights V^^^ are always 
positive; i.e., the sums can be evaluated at arbitrary 
precision, without any sign problem. As in HF-QMC, the 
problem is solved by Monte Carlo importance sampling 
of the auxiliary field {/i} and evaluation of the Green 
function at time slice L with 



G[t^ = [1 + . . . Bi,,Ba,<, . . . B, 



l,<y\ 



(14) 



As a spin flip in the auxiliary field hi(V) at time slice 
I and site i only affects Bi „ at this site, the ratio of 
the weights, needed for the decision whether a proposed 
spin flip is accepted, involves only local quantities; a full 
recomputation of the determinants oi N x N matrices 
appearing in Eq. (|lip is not needed. The computational 
effort is further reduced by calculating the Green function 
at time slice l + l from the quantities at time slice I, using 
so-called "wrapping" : 



Gl+l,a — Bi ^Gl^aBl^cr ■ 



(15) 



In order to avoid the accumulation of numerical errors in 
the matrix multiplications, it is necessary to recalculate 
the full Green function at regular intervals. This is 
particularly important at low temperatures. 

All this considered, the numerical cost scales cubically 
with the number of sites and linearly with the number 



of time slices; at constant Ar, this translates to a total 
effort 0{N^I3), where N = Nb + 1. Note that a need 
for finer bath discretizations at lower temperatures could 
potentially spoil the scaling advantage of the method over 
direct impurity solvers; we will show in Sec. IIII CI that 
this is not the case for our test applications. 

The application in the DMFT context [3l[ starts with 
the computation of the Hamiltonian parameters (for 
some choice of Nt,), exactly like in the ED approach. As 
in the HF-QMC approach, one then chooses some dis- 
cretization At, computes {G(ZAt)}^q for the impurity 
site, and applies a (nontrivial) Fourier transform back to 
Matsubara frequencies. The result is an impurity solver 
with superior scaling (linear in /3) compared to the direct 
impurity solvers (cubic in /3), however with a bias due to 
the Trotter discretization Ar (in addition to a possible 
bias due to the bath discretization with Ni, sites) which, 
as we will show in See. IIIIl can be quite significant. 



C. Specification of multigrid BSS-QMC algorithm 

The central feature of our new algorithm is the elim- 
ination of this systematic Trotter error, while retaining 
the advantage of linear-in-/? scaling inherent in the BSS- 
QMC method. In the following, we will specify the 
method and illustrate it using an example (Fig. [3]), that 
will be discussed in detail in Sec. IIIIl 

In contrast to the previous DMFT-BSS implementa- 
tion with a unique discretization Ar in all BSS com- 
putations throughout the DMFT self-consistency cycle, 
the (impurity) Green function of the Hamiltonian -ffAnd 
at hand is computed in M > 20 parallel BSS runs 
(indexed by 1 < i < M) , each employing a homogeneous 
imaginary-time grid with a specific discretization (AT)i, 
chosen from a set {(AT)i|(Ar)niin < (AT)i < (Ar)i„ax} 
with typically 6-9 different elements. Green functions 
resulting from BSS-QMC runs with the same discretiza- 
tion (At); = (AT)j are averaged over, thereby reducing 
the dependencies on initialization conditions and further 
enhancing the parallelism. 

This leads to a set of Green functions defined, in gen- 
eral, on incommensurate imaginary-time grids (symbols 
in Fig. 13]). In order to apply a local Ar — >• extrapo- 
lation, all G(AT)i have to be transformed to a common 
grid. This is possible since the true G(t) is a smooth 
function; however, a direct spline interpolation of the 
raw QMC results, neglecting higher derivatives, would 
not be accurate [4J]. Instead, we consider differences 
between the raw data {G(Ar)i (^AT)}^g and a reference 
Green function, obtained via Eq. (O from a model self- 
energy EJj?^(iw„) [35l.l!36j, written here for the single-band 
case (for multi-band generalizations, see Ref. ^M,)'- 

K'\ii^n) = U Un^,) - + (n_,) (1 - (n_,)) 



1 



1 



(16) 
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which recovers the exact high-frequency asymptotics of 
I](icj„) and G{iuJn) for any choice of the free parameter 
loq and, therefore, approximates the second and higher- 
order derivatives of G(t) at r — >■ (and t ^ j3) well. This 
match can be further improved by adjusting wq- Con- 
sequently, the differences {G(Ar),(^AT) - G™f(;Ar)}^o 
have smaller absolute values and much smaller higher 
derivatives than the original data; in particular, their 
curvature vanishes asymptotically at the boundaries [i^ . 
Thus, they are well represented by natural cubic splines. 

Usually, the parameters of the piecewise polynomials 
constituting such a spline /spline (2;) are determined from 
discrete data {/mcas(a;i)}i^o such that the discrete data 
are reproduced exactly: /spiine(a;i) = fmeasixi) for all < 
i < N. However, in the QMC context, all measurements 
have statistical errors, i.e., the discrete data are better 
represented as {/meas(xi)±A/ineas(a;i)}iIo '^ith standard 
deviations A/meas, which are also estimated within the 
QMC procedure. It is clear that the usual interpolating 
splines, which do not take the uncertainties of the discrete 
data into account, contain more features than warranted 
by the data (in particular at the Nyquist frequency); in 
the context of Green functions this includes the possi- 
bility of acausal behavior. We use, instead, smoothing 
spline fits [13, El] which reproduce the discrete data 
only within error bars, which are typically 0(10"^), (and 
minimize the curvatures under this constraint); these fits 
can be computed in a very similar procedure and at the 
same cost as interpolating splines. 

After combining these approximating "difference" 
splines with exact expressions for G'^'^^{t) resulting from 
Eqs. ([T5|) and ([S]), we obtain smooth approximations of 
the Green functions, as seen in Fig. [2l^a); the inset also 
demonstrates slight deviations of the continuous spline 
fits from the discrete data (within error bars), e.g., for the 
discretization At = 0.7 (dotted line) at r w 1.4 (circle), 
while most other data points are reproduced within the 
line widths. 

These smooth approximations can be evaluated on an 
arbitrarily fine common grid (e.g. with Arfino = 0.005) 
and extrapolated to At — > 0. This is illustrated in Fig. 
2] for the representative values of r denoted by arrows 
in Fig. Elja). Even though most of the raw BSS-QMC 
data do not include estimates of the Green functions at 
these precise values of r, the transformed data (symbols 
in Fig. 2]) depend very regularly on At, falling on nearly 
straight lines as a function of (Ar)^. Therefore, they 
can accurately and reliably be extrapolated to At — !> 
(lines in Fig. |3] and symbols at At = 0); an application of 
this procedure at all t (on the fine grid) leads to quasi- 
continuous Green functions without significant Trotter 
errors, shown as solid lines in Fig. [3l These results 
can be Fourier transformed to Matsubara frequencies 
in a straightforward manner [cf. Fig. [Ud)]. A similar 
approach has also been useful for computing unbiased 
spectra from BSS-QMC 

At first sight, the computational advantage of the 
multigrid procedure is less obvious in the BSS-QMC 
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FIG. 3. (Color online) BSS-QMC impurity Green functions at 
T = 0.04 (symbols) using a bath representation with A'^;, = 4 
sites (with parameters of converged DMFT-ED solution, long- 
dashed lines) and results of multigrid extrapolation to Ar — 
(solid lines). Upper panel: metallic phase {U = 4.4). Lower 
panel: insulating phase {U = 5.1). Arrows denote r values 
for which the extrapolation is shown in Fig. 21 



context than for HF-QMC [SJ, [S^j since the numerical 
effort for direct computations at small At grows only 
linearly, not cubically, with (At)~^ in the BSS case 
(while the systematic errors decay generically as (At)^ 
for a given impurity problem). However, even for a 
fixed Hamiltonian, so much accuracy can be gained by 
extrapolation that it more than offsets the cost of the 
additional grid points. This is true, in particular, since 
stable results are best obtained by averaging over inde- 
pendent BSS-QMC runs; performing these on variable 
grids then allows for extrapolation without additional 
cost. Furthermore, the individual runs thermalizc faster 
in the multigrid variant, due to the smaller number 
of time slices (and proportionally shorter run time per 
sweep), which enhances the parallelism. Most impor- 
tantly, as we will see below, the DMFT self-consistency 
can magnify any bias of the employed impurity solvers in 
complicated ways (in the vicinity of phase transitions), so 
that controlled results are really dependent on unbiased 
methods, such as our multigrid approach. 
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FIG. 4. (Color online) BSS-QMC estimates of imaginary- 
time Green functions Gat(t) at T = 0.04, U = 4.4 after 
interpolation (corresponding to the colored broken lines in 
Fig. [21 for selected values of r (symbols) and extrapolation to 
At = using least-squares fits (lines). 



III. RESULTS 

In this section, we compare results of the new nu- 
merically exact "multigrid" BSS-QMC method with raw 
BSS-QMC results (at finite Trotter discretization), with 
reference ED results (which are exact at the level of 
the auxiliary Hamiltonian) , and with the predictions of 
established irnpurity solvers (multigrid HF-QMC [13] and 
CT-HYB [2^, [2^, [5l| ) . These comparisons are performed 
in three stages: In Sec. IIII A[ we keep the bath Q and 
its approximation by an auxiliary Hamiltonian fixed and 
discuss the impact of the Trotter error and its elimination 
without the complications of the DMFT self-consistency. 
In Sec. IIII B[ we compare full DMFT solutions obtained 
using the various algorithms at moderate temperature 
(T = 0.04), focusing on the impact of Trotter errors 
on the resulting estimates of double occupancy and 
quasiparticle weight. Finally, we present results also 
at lower temperatures T > 0.01 (with DMFT self- 
consistency) , where the impact of the bath discretization 
becomes particularly relevant, in Sec. IIII Cl 

Following the established practice for the evaluation of 
DMFT impurity solvers (23ll28l|. all of these comparisons 
are performed for the half-filled Hubbard model with 
semi-elliptic "Bethe" density of states [s^] (full band 
width W ~ 4) within the paramagnetic phase. Specifi- 
cally, we choose temperatures T < 0.04, which are below 
the critical temperature T* « 0.055 [la, of the first- 
order metal-insulator transition, and interactions close to 
or within the coexistence region of metallic and insulating 
solutions, which arises from the mean-field character of 
the DMFT. 



A. Green function extrapolation at fixed bath 
Hamiltonian parameters 

In general, a bias present in an impurity solver has a 
two- fold impact: On the one hand, it affects estimates 



of Green functions and all other properties for a given 
impurity problem, defined by its bath Green function Q. 
On the other hand, it shifts the fixed point of the DMFT 
self-consistency cycle, i.e., it also modifies the converged 
bath Green function, which, in turn, also affects the 
measured Green functions and all other properties. In 
this subsection, we study the first effect in isolation 
by fixing the bath Green function to the converged 
solution of the ED procedure for Ni, ~ A bath sites (with 
Hamilton parameters {ii,Vi}f^i). As the same auxiliary 
Hamiltonian is used also in the BSS-QMC algorithm, the 
ED estimates of the Green function are exact for the 
purpose of the current comparison; the impact of the 
bath discretization (which corresponds to a bias on the 
DMFT level) will be discussed in Sec. IhTCI 

Local imaginary-time Green functions G(t) are shown 
in Fig. ^a.) for the metallic phase, at J7 = 4.4, and in 
Fig. |3Jb) for the insulating phase, at C/ = 5.1. Here and 
in the following, we restrict the imaginary-time range 
to < T < 13/2; data for t > 13/2 follow from the 
particle- hole symmetry G'(/3 — t) = G(r). Symbols (in 
the magnified insets) represent raw BSS-QMC results 
(with discretizations At = 0.4, At = 0.7, and At = 
0.9); colored long-dashed, dotted, and dash-dotted lines 
denote interpolations obtained using the methods de- 
scribed in Sec. Ill Cl Due to the large discretization, these 
data deviate significantly from the ED reference results 
(gray long-dashed lines), in particular at moderately low 
imaginary times t w 2. In contrast, multigrid BSS-QMC 
Green functions (solid lines) are indistinguishable from 
the ED data at [/ = 5.1 and very close to them at 
U = 4.4, with deviations of the order of statistical errors. 
Thus, our method yields, indeed, quasi-continuous Green 
functions without significant Trotter errors in both test 
cases, although the discretizations of the underlying raw 
BSS-QMC computations (with 0.3 < At < 1.0) would be 
considered much too coarse in conventional applications. 

A very similar picture emerges in an analogous compar- 
ison for the two coexisting solutions aX U = 4.74, shown 
in Fig. [5j Again, the raw BSS-QMC results (symbols 
and colored broken lines) show a strong systematic bias, 
towards more metallic Green functions and of different 
magnitude in the different phases, while the extrapolated 
Green functions agree nearly perfectly with the ED 
references. In fact, some of the BSS Green functions 
calculated for an insulating bath (lower set of symbols 
and broken lines) show such large discretization errors 
at small t < 2, that they approach the exact Green 
function of the metallic DMFT solution (upper solid 
and long-dashed lines). One may suspect from this 
observation that these biased "insulating" solutions will 
not be associated with stable DMFT fixed points if they 
are fed back in the self-consistency cycle; such shifts of 
stability regions induced by the Trotter bias at At > 
will, indeed, be seen in Sec. IIII Bl 

It is clear that the proposed multigrid extrapolation 
technique can only be useful as a practical method if it is 
insensitive to the particular set {(At)^} of discretizations 
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FIG. 5. (Color online) BSS-QMC impurity Green functions 
at T = 0.04 and U = 4.74 (symbols and colored broken 
lines) using a bath representation with Nt = 4 sites (with 
parameters fixed by converged DMFT-ED solution, long- 
dashed lines) and extrapolation to At — (solid lines). 
Upper (lower) set of curves: metallic (insulating) bath. 
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FIG. 6. (Color online) Green functions in the insulating 
phase at T = 0.04, U = 4.7, extrapolated from BSS-QMC 
results using different imaginary-time grids [53] (at fixed bath 
representation). The excellent agreement shows that the 
multigrid procedure is stable with respect to its technical 
parameters. 

in the underlying BSS-QMC runs, i.e., if no sensible 
choice leads to a significant bias. This is demonstrated 
in Fig. [S] for the insulating phase at U = 4.7: the Green 
functions for the same auxiliary problem obtained from 
multigrid extrapolations with three different At grids 
[53! agree perfectly within the precision of the method. 
The latter is primarily determined by the statistical 



errors, i.e., by the number of sweeps and, possibly, 
by the numerical precision in the matrix operations. 
Only if raw BSS-QMC data of much higher precision 
was available (with many millions of sweeps per run), 
additional accuracy could be gained by choosing smaller 
discretizations (e.g., 0.1 < At < 0.3). As a rule 
of thumb, the multigrid procedure can be based on 
discretizations At that are 3 to 10 times as large as the 
discretization that one would choose in a conventional 
BSS-QMC procedure. 



B. Comparisons of impurity solvers at full DMFT 
self-consistency: impact of Trotter errors 

So far, we have compared different algorithms just 
at the impurity level, i.e., for a fixed bath Green 
function (determined from a self-consistent DMFT-ED 
calculation). In contrast, we will now discuss results of 
completely independent DMFT solutions, each of which 
corresponds to full self-consistency for a given impurity 
solver (cf. Fig. [5]). For all Hamiltonian based methods 
(ED, BSS-QMC, and muhigrid BSS-QMC), the number 
of bath sites is restricted to A^f, = 4 (as above) ; the impact 
of this parameter will be studied in Sec. IIIICI 

Specifically, we discuss static observables that are 
particularly useful for discriminating between metallic 
and (possibly coexisting) insulating DMFT solutions, 
namely the double occupancy 



D = {n^ni) , 



(17) 



which is proportional to the interaction energy i?int 
UD, and the quasiparticle weight 





aReE(w) 




-1 


z = 






did 


U1 = 





Im E(2Wi) 



1 -1 



ttT 



(18) 

Open symbols in Fig.[7]denote estimates resulting from 
self-consistent DMFT solutions using the conventional 
BSS-QMC impurity solver at finite discretization 0.3 < 
At < 0.5, i.e., using the scheme established in Ref. [sj 
The estimated values of Z, shown in Fig. [Tljb), have a 
nearly uniform offset in the metallic phase at U < 4.8 
relative to each other and relative to the reference ED 
solution (gray diamonds). The Trotter bias inherent 
in the conventional BSS-QMC procedure also leads to 
a significant overestimation of the range of stability of 
the metallic solution: The metallic BSS-QMC solutions 
extend to much larger interactions (e.g., to U ^ 5.1 at 
At = 0.4) than the ED reference solution. 

This is also seen in corresponding estimates of the 
double occupancy [Fig. El^a)]; however, for these ob- 
servables the Trotter bias is highly nonuniform (in the 
metallic solution): a.t U = 4.7 (arrow), the conventional 
BSS-QMC estimates are nearly on top of each other; 
relative deviations are only clearly seen at stronger 
interactions U > 4.9 and (to a lesser degree) at weaker 
interactions U < 4.5. At the same time, nearly all 
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FIG. 8. (Color online) Discretization dependence of the 
double occupancy D as estimated from BSS-QMC, using 
either the multigrid scheme (filled symbols) or self-consistent 
BSS-QMC solutions at finite At (open symbols), within the 
metallic phase at f/ = 4.7 (upper data set) or the insulating 
phase at !7 = 5.0 (lower data set). Dashed (solid) lines denote 
least squares fits to the multigrid (conventional) BSS-QMC 
data at At > 0.3; dotted lines denote fits that include also 
data at At = 0.1 and At = 0.2 (small open symbols). 



FIG. 7. (Color online) Estimates of double occupancy D{U) 
and quasiparticle weight Z{U) obtained in independent self- 
consistent DMFT calculations using various impurity solvers: 
multigrid HF-QMC (crosses), conventional BSS-QMC (open 
symbols), multigrid BSS-QMC (circles), and ED (diamonds). 
In each panel, the upper (lower) sets of curves correspond 
to metallic (insulating) solutions. Lines are guides to the 
eye only. Arrows in (a) indicate parameters for which the 
discretization dependence is studied in Fig. [S] 

of these data deviate significantly (and without obvious 
systematics) from the reference ED result (diamonds) , so 
that an a posteriori elimination of the Trotter bias seems 
impossible. 

In contrast, the new multigrid BSS-QMC procedure, as 
discussed in Sec. Ill Cl and illustrated in Fig.[21Jd), leads to 
estimates of both D and Z (filled circles) which perfectly 
recover the ED solutions, even though they arc based on 
BSS-QMC runs with At > 0.3. 

This is also true for the insulating solutions (lower 
sets of curves in Fig. [7]), the stability range of which 
is also shifted towards stronger interactions in the case 
of conventional BSS-QMC calculations (open symbols); 
here, the Trotter bias appears roughly uniform for D and 
very nonuniform for Z . Again, the multigrid BSS-QMC 
results agree perfectly with the ED reference data. 

For comparison, crosses and black solid lines in Fig. [7] 
denote estimates of an unbiased direct impurity solver, 
namely the multigrid HF-QMC method (sj; these show 
good overall agreement with both the ED and the 
multigrid BSS-QMC data. A slight negative deviation 
in the estimates of D of the latter, Hamiltonian based, 
methods can be traced back to the relatively poor bath 
discretization with iV^ = 4 auxiliary sites (cf. Sec. IIII C|) . 

Since the double occupancy D is best computed di- 
rectly on the impurity level (in QMC based approaches) , 



its physical value has to be extrapolated from raw 
estimates I?Ar, with discretizations corresponding to the 
different grid points used within the multigrid procedure 
(in contrast to the quasiparticle weight Z, which follows 
from the self-energy E, which, in turn, is determined from 
unbiased Green functions). 

As seen in Fig. [51 the Trotter bias inherent in these raw 
estimates (filled symbols) is perfectly regular Q even at 
large At, so that reliable extrapolations At — > (thick 
dashed lines) are possible both in the metallic phase, 
at [/ = 4.7 (upper set of curves), and in the insulating 
phase, at [/ ~ 5.0 (lower set of curves). 

In contrast, estimates oi D resulting from conventional 
BSS-QMC calculations in the same range of discretiza- 
tions At > 0.3 (large open symbols in Fig. ^ show 
such irregular dependencies on At that quadratic least- 
square fits (solid lines) lead to extrapolations At 
with significant offsets. Roughly accurate results (dotted 
lines) can only be obtained when including raw data at 
much smaller discretizations (small open symbols). This 
shows, again, that only an elimination of all Trotter 
errors within the self-consistency cycle, as introduced 
by our multigrid approach, can efficiently generate high- 
precision results. 

C. Comparisons of impurity solvers at full DMFT 
self-consistency: impact of bath discretization 

So far, we have restricted the bath representation in 
all Hamiltonian based impurity solvers (ED and both 
variants of BSS-QMC) to only iVj, = 4 bath sites and 
focused on the impact of the Trotter errors and their 
elimination. From the mutual agreement with multigrid 
HF-QMC, an impurity solver which treats the bath 
directly on the action level, we can conclude that this 
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FIG. 9. (Color online) Estimates of double occupancy D{U) 
and quasiparticle weight Z{U) at T = 0.04, obtained in 
selfconsistent DMFT calculations using Hamiltonian based 
impurity solvers with 3 < A^;, < 6 bath sites: multigrid 
BSS-QMC (chcles), ED (open symbols). Multigrid HF-QMC 
results (crosses) represent the limit Nt ^ ca. In each 
panel, the upper (lower) sets of curves correspond to metallic 
(insulating) solutions. Lines are guides to the eye only. 



coarse bath discretization allows for reasonably accurate 
estimates of D and, in particular, Z at the moderately 
low temperature T = 0.04. However, the ED and 
multigrid BSS-QMC estimates of D were found in Fig. 
[7] to lie a bit below the multigrid HF-QMC data; this 
deviation must be an artifact of the bath discretization 
if the multigrid HF-QMC reference data are correct. 
Moreover, we must suspect that the bath discretization 
bias gets worse (at constant Ni,) at lower temperatures. 

Figure IH] shows estimates of D{U) and Z{U) at T = 
0.04, similarly to Fig. [7] and with the same multigrid 
HF-QMC reference data (crosses), but now using Hamil- 
tonian based impurity solvers with 3 < A'^t, < 6 bath 
sites. Here and in the following, "BSS" refers to multigrid 
BSS-QMC data, i.e., without significant Trotter errors; 
for simplicity, we have used this method only for the 
largest auxiliary Hamiltonian {Nh ~ 6). Smaller bath 
sizes {Nb = 3, Nb = 4, and = 5) are represented 
only by the ED solution, which is cheaper and free of 
statistical noise. At the resolution of Fig.[9l the estimates 
associated with the finer bath discretizations Ni, = 5 
(triangles) and TV;, = 6 (circles) agree with each other. 
Therefore and since they are also consistent with the 
unbiased multigrid HF-QMC data (crosses), we conclude 
that convergence with respect to the bath discretization 
is reached already at A'b = 5 at T = 0.04. In contrast, 
the ED estimates of D are apparently slightly too small 
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FIG. 10. (Color online) Estimates of double occupancy D{U) 
and quasiparticle weight Z{U) at T = 0.02, using bath 
discretizations with 4 < A*'(, < 6 sites, analogous to Fig. [9] 



at TVb = 4 (pentagons); corresponding results at A^t = 3 
(squares) arc far off both for D and Z . 

Note that consistent convergence of observable esti- 
mates with Nb; as demonstrated in Fig. [9] (as well as 
Fig. [TU] and Fig. [TT|) can only be observed when optimal 
Hamiltonian parameters are determined with great care, 
as described in Sec. Ill A[ within each self-consistency 
cycle; otherwise some bath sites may remain ineffective or 
the estimates can even get worse upon increasing Nb- In 
addition (as always in the DMFT context), it is essential 
that enough DMFT iterations are performed at each 
phase point in order to ensure convergency with respect 
to the self-consistency cycle (cf. Fig. [5]). 

Halving the temperature amplifies the bath discretiza- 
tion effects, as seen in Fig. (TU] At T = 0.02, only 
the best Hamiltonian representation [Nb — 6, evaluated 
with multigrid BSS-QMC, circles) recovers all reference 
multigrid HF-QMC results (crosses) within their accu- 
racy. At TVb = 5, the estimates of D in the insulating 
phase are already slightly too small; at Nb — 4, strong 
negative deviations in D{U) are apparent also for the 
metallic solution. The impact of the bath discretization 
becomes even much stronger at T = 0.01, as shown in 
Fig. [TT] for the metallic phase, which is interesting as a 
strongly renormalized Fermi liquid (while the properties 
of the insulating phase are asymptotically independent of 
temperature). We find that even the results for Nb = 5 
and Nb = 6 deviate significantly in Fig. ITlT a) from each 
other and from the multigrid HF-QMC reference result, 
especially at C/ = 5.3, near the edge of the stability 
region of the metallic phase. Only the multigrid BSS- 
QMC results using Nb = 7 bath sites (circles) agree 
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FIG. 11. (Color online) Estimates of double occupancy D{U) 
and quasiparticle weight Z{U) at T = 0.01, using bath 
discretizations with 4 < A^i, < 7 sites, analogous to Fig. 
[51 The CT-HYB data (diamonds) also represent the limit 
Nb oo. 

with the reference data within their precision; even better 
agreement is observed with data obtained using the CT- 
HYB impurity solver (diamonds). Note that BSS-QMC 
is much more efficient than ED at Ni, = 7; in particular, 
the latter would need orders of magnitude more main 
memory than the former. 

As noted in Sec. Ill Bl a strong increase with inverse 
temperature of the number Nb of bath sites needed 
for a given accuracy could, in principle, eliminate the 
scaling advantage of the DMFT-BSS approach, as its 
computational cost does not only include the direct 
factor but also a factor = N^{/3). However, 
our results indicate that this effect is minor: In our 
test case, we needed to add one bath site upon halving 
the temperature for roughly constant accuracy; this is 
consistent with a scaling Nb cx ln(/3), i.e., an overall 
computational cost proportional to /3[ln(/3)]^ which is 
still linear up to logarithmic corrections. 



IV. CONCLUSIONS 

The DMFT and its extensions are invaluable tools for 
the study of phenomena associated with strong electronic 
correlations and for quantitative predictions of properties 
of correlated materials. However, the numerical solution 
of the DMFT self-consistency equations remains a great 
challenge: the established, direct, QMC impurity solvers 
yield unbiased results, but provide only limited access 
to the low-T phase regions of interest, due to the cubic 
scaling of their computational cost with the inverse 
temperature /?. Exact diagonalization (ED) approaches, 
on the other hand, are limited by their exponential 
scaling with the number of sites N of the auxiliary 
Hamiltonian. 

The multigrid BSS-QMC algorithm presented in this 
work allows for solving the DMFT self-consistency equa- 
tions with an effort that grows only linearly with /3; 
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in contrast to an earlier BSS-QMC based method 
it is free of significant Trotter errors, i.e., numerically 
exact at the level of the auxiliary Hamiltonian. Since the 
computational cost grows only cubically with N, much 
better representations of the bath are possible than for 
ED. As demonstrated by applications to the half- filled 
Hubbard model in and near the coexistence region of 
metallic and insulating solutions and by comparisons 
with direct QMC impurity solvers, the new method yields 
unbiased results (for sufficiently ffiie bath discretization), 
in spite of using quite coarse Trotter discretizations in the 
underlying BSS-QMC evaluations. 

The new unbiased quasi-CT impurity solver should 
show its full potential in multi-band cases and in cluster 
extensions of DMFT, where the prefactor A^'^ of the BSS- 
QMC scheme (compared to a factor of 1 in HF-QMC 
calculations in single-site DMFT) is leveled off by the 
increased complexity of the original DMFT problem. Our 
approach can also be extended beyond Hubbard models; 
it could be particularly valuable for the cellular DMFT 
treatment of the Kondo lattice model, where interesting 
temperature regimes are out of reach of the existing 
impurity solvers j26l [ssj . 
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